Analiza zbioru była dosyć skomplikowana z powodu dużej liczby kolumn.
Dodatkowo nie ułatwiała tego dość egzotyczna, biologiczna dziedzina.
Problemy pojawiły się przy określaniu korelacji pomiędzy kolumnami, gdzie część nie była numeryczna, a poza tym, nie sposób było pokazać wszystkich równocześnie.
library(dplyr)
library(tibble)
library(tidyr) #transformowanie danych
library(ggplot2)
library(plotly)
library(kableExtra) #ładne tabele
library(data.table) #%like%
library(corrplot)
library(caret)
set.seed(23)
Wczytujemy dane z pliku all_summary.csv. Za pomocą heurystyki wyznaczamy klasy pierwszych 100 wierszy.
filePath = "C:/Users/micha/Downloads/all_summary.csv"
initial <- read.table(filePath, nrows = 100, sep = ";", header = TRUE)
classes <- sapply(initial, class)
file_data <- read.table(filePath, colClasses = classes, sep = ";", na.strings = 'nan', nrows = 10000, header = TRUE)
data <- file_data
Dane znajdują się w zmiennej data.
Liczba wierszy w zbiorze: 10000.
Usuwamy wiersze posiadające wartość zmiennej ref_name równą: “UNK”, “UNX”, “UNL”, “DUM”, “N”, “BLOB”, “ALA”, “ARG”, “ASN”, “ASP”, “CYS”, “GLN”, “GLU”, “GLY”, “HIS”, “ILE”, “LEU”, “LYS”, “MET”, “MSE”, “PHE”, “PRO”, “SEC”, “SER”, “THR”, “TRP”, “TYR”, “VAL”, “DA”, “DG”, “DT”, “DC”, “DU”, “A”, “G”, “T”, “C”, “U”, “HOH”, “H20”, “WAT”
res_names_to_remove <- c('UNK', 'UNX', 'UNL', 'DUM', 'N', 'BLOB', 'ALA', 'ARG', 'ASN', 'ASP', 'CYS', 'GLN', 'GLU', 'GLY', 'HIS', 'ILE', 'LEU', 'LYS', 'MET', 'MSE', 'PHE', 'PRO', 'SEC', 'SER', 'THR', 'TRP', 'TYR', 'VAL', 'DA', 'DG', 'DT', 'DC', 'DU', 'A', 'G', 'T', 'C', 'U', 'HOH', 'H20', 'WAT')
data <- filter(data, !res_name %in% res_names_to_remove)
Liczba wierszy w zbiorze po usunięciu wierszy z ww. wartościami res_name: 9940.
W poniższej sekcji oczyścimy dane z kolumn i wierszy zawierających wartości NaN.
Na początku analiza kolumn.
Wyznaczamy kolumny, które zawierają przynajmniej jedną wartość NaN, grupujemy po nazwie.
column_value_DT <- gather(data, "column_name", "val", 1:ncol(data))
columns_with_nans <- filter(column_value_DT, is.na(column_value_DT$val))[,1]
Wyświetlamy podsumowanie:
table(columns_with_nans)
## columns_with_nans
## dict_atom_C_count dict_atom_N_count
## 174 174
## dict_atom_non_h_count dict_atom_non_h_electron_sum
## 174 174
## dict_atom_O_count dict_atom_S_count
## 174 174
## part_01_density_CI part_01_density_E2_E1
## 80 80
## part_01_density_E3_E1 part_01_density_E3_E2
## 80 80
## part_01_density_FL part_01_density_FL_norm
## 80 80
## part_01_density_I1 part_01_density_I1_norm
## 80 80
## part_01_density_I2 part_01_density_I2_norm
## 80 80
## part_01_density_I3 part_01_density_I3_norm
## 80 80
## part_01_density_I4 part_01_density_I4_norm
## 80 80
## part_01_density_I5 part_01_density_I5_norm
## 80 80
## part_01_density_I6 part_01_density_I6_norm
## 80 80
## part_01_density_M000 part_01_density_O3
## 80 80
## part_01_density_O3_norm part_01_density_O4
## 80 80
## part_01_density_O4_norm part_01_density_O5
## 80 80
## part_01_density_O5_norm part_01_density_sqrt_E1
## 80 80
## part_01_density_sqrt_E2 part_01_density_sqrt_E3
## 80 80
## part_01_density_Z_0_0 part_01_density_Z_1_0
## 80 80
## part_01_density_Z_2_0 part_01_density_Z_2_1
## 80 80
## part_01_density_Z_3_0 part_01_density_Z_3_1
## 80 80
## part_01_density_Z_4_0 part_01_density_Z_4_1
## 80 80
## part_01_density_Z_4_2 part_01_density_Z_5_0
## 80 80
## part_01_density_Z_5_1 part_01_density_Z_5_2
## 80 80
## part_01_density_Z_6_0 part_01_density_Z_6_1
## 80 80
## part_01_density_Z_6_2 part_01_density_Z_6_3
## 80 80
## part_01_density_Z_7_0 part_01_density_Z_7_1
## 80 80
## part_01_density_Z_7_2 part_01_density_Z_7_3
## 80 80
## part_01_shape_CI part_01_shape_E2_E1
## 80 80
## part_01_shape_E3_E1 part_01_shape_E3_E2
## 80 80
## part_01_shape_FL part_01_shape_FL_norm
## 80 80
## part_01_shape_I1 part_01_shape_I1_norm
## 80 80
## part_01_shape_I2 part_01_shape_I2_norm
## 80 80
## part_01_shape_I3 part_01_shape_I3_norm
## 80 80
## part_01_shape_I4 part_01_shape_I4_norm
## 80 80
## part_01_shape_I5 part_01_shape_I5_norm
## 80 80
## part_01_shape_I6 part_01_shape_I6_norm
## 80 80
## part_01_shape_M000 part_01_shape_O3
## 80 80
## part_01_shape_O3_norm part_01_shape_O4
## 80 80
## part_01_shape_O4_norm part_01_shape_O5
## 80 80
## part_01_shape_O5_norm part_01_shape_sqrt_E1
## 80 80
## part_01_shape_sqrt_E2 part_01_shape_sqrt_E3
## 80 80
## part_01_shape_Z_0_0 part_01_shape_Z_1_0
## 80 80
## part_01_shape_Z_2_0 part_01_shape_Z_2_1
## 80 80
## part_01_shape_Z_3_0 part_01_shape_Z_3_1
## 80 80
## part_01_shape_Z_4_0 part_01_shape_Z_4_1
## 80 80
## part_01_shape_Z_4_2 part_01_shape_Z_5_0
## 80 80
## part_01_shape_Z_5_1 part_01_shape_Z_5_2
## 80 80
## part_01_shape_Z_6_0 part_01_shape_Z_6_1
## 80 80
## part_01_shape_Z_6_2 part_01_shape_Z_6_3
## 80 80
## part_01_shape_Z_7_0 part_01_shape_Z_7_1
## 80 80
## part_01_shape_Z_7_2 part_01_shape_Z_7_3
## 80 80
## part_02_density_CI part_02_density_E2_E1
## 704 704
## part_02_density_E3_E1 part_02_density_E3_E2
## 704 704
## part_02_density_FL part_02_density_FL_norm
## 704 704
## part_02_density_I1 part_02_density_I1_norm
## 704 704
## part_02_density_I2 part_02_density_I2_norm
## 704 704
## part_02_density_I3 part_02_density_I3_norm
## 704 704
## part_02_density_I4 part_02_density_I4_norm
## 704 704
## part_02_density_I5 part_02_density_I5_norm
## 704 704
## part_02_density_I6 part_02_density_I6_norm
## 704 704
## part_02_density_M000 part_02_density_O3
## 704 704
## part_02_density_O3_norm part_02_density_O4
## 704 704
## part_02_density_O4_norm part_02_density_O5
## 704 704
## part_02_density_O5_norm part_02_density_sqrt_E1
## 704 704
## part_02_density_sqrt_E2 part_02_density_sqrt_E3
## 704 704
## part_02_density_Z_0_0 part_02_density_Z_1_0
## 704 704
## part_02_density_Z_2_0 part_02_density_Z_2_1
## 704 704
## part_02_density_Z_3_0 part_02_density_Z_3_1
## 704 704
## part_02_density_Z_4_0 part_02_density_Z_4_1
## 704 704
## part_02_density_Z_4_2 part_02_density_Z_5_0
## 704 704
## part_02_density_Z_5_1 part_02_density_Z_5_2
## 704 704
## part_02_density_Z_6_0 part_02_density_Z_6_1
## 704 704
## part_02_density_Z_6_2 part_02_density_Z_6_3
## 704 704
## part_02_density_Z_7_0 part_02_density_Z_7_1
## 704 704
## part_02_density_Z_7_2 part_02_density_Z_7_3
## 704 704
## part_02_shape_CI part_02_shape_E2_E1
## 704 704
## part_02_shape_E3_E1 part_02_shape_E3_E2
## 704 704
## part_02_shape_FL part_02_shape_FL_norm
## 704 704
## part_02_shape_I1 part_02_shape_I1_norm
## 704 704
## part_02_shape_I2 part_02_shape_I2_norm
## 704 704
## part_02_shape_I3 part_02_shape_I3_norm
## 704 704
## part_02_shape_I4 part_02_shape_I4_norm
## 704 704
## part_02_shape_I5 part_02_shape_I5_norm
## 704 704
## part_02_shape_I6 part_02_shape_I6_norm
## 704 704
## part_02_shape_M000 part_02_shape_O3
## 704 704
## part_02_shape_O3_norm part_02_shape_O4
## 704 704
## part_02_shape_O4_norm part_02_shape_O5
## 704 704
## part_02_shape_O5_norm part_02_shape_sqrt_E1
## 704 704
## part_02_shape_sqrt_E2 part_02_shape_sqrt_E3
## 704 704
## part_02_shape_Z_0_0 part_02_shape_Z_1_0
## 704 704
## part_02_shape_Z_2_0 part_02_shape_Z_2_1
## 704 704
## part_02_shape_Z_3_0 part_02_shape_Z_3_1
## 704 704
## part_02_shape_Z_4_0 part_02_shape_Z_4_1
## 704 704
## part_02_shape_Z_4_2 part_02_shape_Z_5_0
## 704 704
## part_02_shape_Z_5_1 part_02_shape_Z_5_2
## 704 704
## part_02_shape_Z_6_0 part_02_shape_Z_6_1
## 704 704
## part_02_shape_Z_6_2 part_02_shape_Z_6_3
## 704 704
## part_02_shape_Z_7_0 part_02_shape_Z_7_1
## 704 704
## part_02_shape_Z_7_2 part_02_shape_Z_7_3
## 704 704
## weight_col
## 9940
Jak widać na podsumowaniu powyżej, dla wszystkich wierszy kolumna weight_col ma wartość NaN. Usuwamy ją ze zbioru danych:
data <- data[,-which(names(data) == 'weight_col' )]
Widać również, że część wierszy nie posiada wartości dla kolumn zaczynających się od part_02 oraz part_01.
W celu oczyszczenia danych usuwamy te wiersze.
all_rows_count <- nrow(data)
data <- data[complete.cases(data),]
cleaned_rows_count <- nrow(data)
removed_rows_count <- all_rows_count - cleaned_rows_count
Usunięto 869 wierszy. Stanowi to 8.7424547 % wszystkich wierszy.
Liczba wierszy: 9071
Liczba kolumn: 411
Liczba unikalnych wartości res_name: 801
Poniżej rozmiar data frame w pamięci (w Mb):
print(object.size(data), units='Mb')
## 36.8 Mb
W tej sekcji filtrujemy dane do wierszy, których klasa znajduje się wśród 50 najpopularniejszych.
Na początku grupujemy dane po atrybucie res_name i sortujemy.
top50 <- sort(table(data$res_name), decreasing=TRUE)[1:50]
top50names <- names(top50)
Poniżej znajduje się posortowana lista 50 najpopularniejszych klas.
top50names
## [1] "SO4" "GOL" "EDO" "NAG" "CL" "DMS" "ZN" "CA" "HEM" "MG" "PO4"
## [12] "NA" "NAD" "ACT" "FAD" "PEG" "IOD" "CD" "K" "PG4" "MN" "NAP"
## [23] "FMN" "MLY" "LHG" "AMP" "NDP" "FE2" "EPE" "ADP" "PLP" "ACY" "HEC"
## [34] "BR" "MAN" "CU" "FE" "FMT" "MPD" "C8E" "CIT" "CME" "PGE" "FES"
## [45] "H4B" "PLC" "CYC" "FEC" "GLC" "NI"
Filtrujemy za pomocą wyznaczonych 50 nazw
data <- filter(data, res_name %in% top50names)
data$res_name <- droplevels(data$res_name)
Pozostało 6376 wierszy.
Korelację wyznaczamy dla kolumn numerycznych.
Pomijamy kolumny: “blob_coverage” “res_coverage” “title” “pdb_code” “res_name” “chain_id” “skeleton_data” “fo_col” “fc_col”
Dla zbioru danych korelację wyznaczamy za pomocą funkcji cor.
Następnie prezentujemy wyniki, za pomocą funkcji corrplot.
is_numeric_column_list <- unlist(lapply(data, is.numeric))
numeric_data <- data[, is_numeric_column_list]
cor_data <- cor(numeric_data, use = "pairwise.complete.obs")
#pdf(file = "korelacja_miedzy_kolumnami_FULL.pdf")
corrplot(cor_data, method = "color", tl.col="black",
tl.cex = 0.2, cl.cex = 0.2, title = "Korelacja między kolumnami")
#dev.off()
Z racji dużego rozmiaru macierzy (402 x 402) wykres jest nieczytelny.
Wynik można zobaczyć również w wyeksportowanym pliku pdf: korelacja_miedzy_kolumnami.pdf
W celu zwiększenia widoczności, ogarniczymy liczbę kolumn, tylko do tych zaczynających się od part_01:
is_part_01_column <- colnames(numeric_data) %like% 'part_01'
part_01_column_names <- colnames(numeric_data)[is_part_01_column]
filtered_numeric_data <- numeric_data[is_part_01_column]
cor_data2 <- cor(filtered_numeric_data, use = "pairwise.complete.obs")
#pdf(file = "korelacja_miedzy_kolumnami_part_01.pdf")
corrplot(cor_data2, tl.col="black", order = "FPC",
tl.cex = 0.3, title = "Korelacja między kolumnami")
#dev.off()
W pliku “korelacja_miedzy_kolumnami_part_01.pdf” znajduje się powyższa wizualizacja.
Dodatkowo sprawdzamy korelację dla kolumn nie zaczynających się od part (których jest najwięcej), oraz skeleton.
no_part_skeleton_data <- numeric_data %>% select(-starts_with("part")) %>% select(-starts_with("skeleton"))
cor_data3 <- cor(no_part_skeleton_data, use = "pairwise.complete.obs")
cor_data3[is.na(cor_data3)] <- 0
#pdf(file = "korelacja_miedzy_kolumnami_no_part_skeleton.pdf")
corrplot(cor_data3, tl.col="black", order = "FPC", method = "square",
tl.cex = 0.5, title = "Korelacja między kolumnami")
#dev.off()
Tutaj z racji dużo mniejszej liczby kolumn, wykres jest bardziej czytelny.
Jak można zauważyć, większość kolumn zawierających w swojej nazwie słowo atom, jest silnie skorelowana.
sort(table(data$res_name), decreasing=TRUE)
##
## SO4 GOL EDO NAG CL DMS ZN CA HEM MG PO4 NA NAD ACT FAD PEG IOD CD
## 982 587 466 396 380 327 318 282 234 194 169 161 129 116 112 104 96 87
## K PG4 MN NAP FMN MLY LHG AMP NDP FE2 EPE ADP PLP ACY HEC BR MAN CU
## 76 71 59 55 54 54 52 45 45 43 42 41 41 37 36 33 33 32
## FE FMT MPD C8E CIT CME PGE FES H4B PLC CYC FEC GLC NI
## 32 31 29 28 28 28 28 27 27 27 26 26 25 25
chart_data_1 <- ggplot(data, aes(x = data$local_res_atom_non_h_count))
chart_layout_1 <- geom_histogram(col = 'blue', fill = 'orange', binwidth = 1)
chart1 <- chart_data_1 + chart_layout_1 + labs(title="Rozkład liczby atomów") + labs(x="Liczba atomów")
chart1
chart_data_2 <- ggplot(data, aes(x = data$local_res_atom_non_h_electron_sum))
#chart_layout_2 <- geom_histogram(col = 'blue', fill = 'orange', binwidth = 5)
chart2 <- chart_data_2 + chart_layout_1 + labs(title="Rozkład liczby elektonów") + labs(x="Liczba elektronów")
chart2
Na początek wybieramy interesujące nas kolumny do porównań: local_res_atom_non_h_count, oraz dict_atom_non_h_count. Dodatkowo każdy wiersz zawiera nazwę klasy.
atoms_incompatibility <- unique(select(data, res_name, local_res_atom_non_h_count, dict_atom_non_h_count))
Dodajemy kolumnę z różnicą, wyliczoną na postawie 2 innych.
atoms_incompatibility_with_diff <- mutate(atoms_incompatibility, diff = abs(local_res_atom_non_h_count - dict_atom_non_h_count))
Sortujemy i wyświetlamy 5 pierwszych.
atoms_incompatibility_with_diff_sorted <- atoms_incompatibility_with_diff[order(-atoms_incompatibility_with_diff$diff),]
head(atoms_incompatibility_with_diff_sorted, 5)
## res_name local_res_atom_non_h_count dict_atom_non_h_count diff
## 63 LHG 8 49 41
## 62 LHG 11 49 38
## 77 PLC 7 42 35
## 78 PLC 9 42 33
## 4 NDP 22 48 26
Poniżej znajduje się tabela pokazująca 10 klas z największą niezgodnością liczby atomów:
atoms_table_data <- select(atoms_incompatibility_with_diff_sorted, res_name, diff)[1:10,]
row.names(atoms_table_data) <- NULL
atoms_table_data %>% knitr::kable(col.names = c('res_name', 'Niezgodność liczy atomów')) %>% kable_styling()
| res_name | Niezgodność liczy atomów |
|---|---|
| LHG | 41 |
| LHG | 38 |
| PLC | 35 |
| PLC | 33 |
| NDP | 26 |
| PLC | 24 |
| PLC | 23 |
| NAP | 22 |
| PLC | 22 |
| NAP | 17 |
Analogicznie postępujemy aby pokazać niezgodność liczby elektronów:
electrons_incompatibility <- unique(select(data, res_name, local_res_atom_non_h_electron_sum, dict_atom_non_h_electron_sum))
electrons_incompatibility_with_diff <- mutate(electrons_incompatibility, diff = abs(local_res_atom_non_h_electron_sum - dict_atom_non_h_electron_sum))
electrons_incompatibility_with_diff_sorted <- electrons_incompatibility_with_diff[order(-electrons_incompatibility_with_diff$diff),]
electrons_table_data <- select(electrons_incompatibility_with_diff_sorted, res_name, diff)[1:10,]
row.names(electrons_table_data) <- NULL
electrons_table_data %>% knitr::kable(col.names = c('res_name', 'Niezgodność liczy elektronów')) %>% kable_styling()
| res_name | Niezgodność liczy elektronów |
|---|---|
| LHG | 275 |
| LHG | 257 |
| PLC | 236 |
| PLC | 224 |
| NDP | 170 |
| PLC | 170 |
| PLC | 164 |
| NAP | 158 |
| PLC | 158 |
| NAP | 112 |
Na początek wyszukujemy nazwy kolumn zaczynające się od part_01. Wyświetlamy 10 pierwszych.
is_part_01_column <- colnames(data) %like% 'part_01'
part_01_column_names <- colnames(data)[is_part_01_column]
head(part_01_column_names, 10)
## [1] "part_01_shape_segments_count" "part_01_density_segments_count"
## [3] "part_01_volume" "part_01_electrons"
## [5] "part_01_mean" "part_01_std"
## [7] "part_01_max" "part_01_max_over_std"
## [9] "part_01_skewness" "part_01_parts"
Następnie filtrujemy dane, przekształcamy do postaci:
filtered_data <- data[is_part_01_column]
gathered_data <- gather(filtered_data, 'col_name', "value", part_01_column_names)
dim(gathered_data)
## [1] 675856 2
wyświetlamy rozkład wartości, osobny wykres dla każdej kolumny:
k <- ggplot(gathered_data, aes(x = gathered_data$value))
i <- k + facet_wrap(vars(gathered_data$col_name), ncol = 3, scales = "free")
i + geom_histogram(bins = 25) + xlab("Rozkład wartości")
p <- ggplot(data, aes(local_res_atom_non_h_count, local_res_atom_non_h_electron_sum,
color=res_name)) +
geom_point() +
xlab("Liczba atomów") +
ylab("Liczba elektonów")
ggplotly(p)
Filtrujemy dane do kolumn numerycznych. Wykorzystujemy funkcję lm.
is_numeric_column_list <- unlist(lapply(data, is.numeric))
numeric_data <- data[, is_numeric_column_list]
without_electron <- names(numeric_data)[-which(names(numeric_data) =='local_res_atom_non_h_electron_sum')]
joined_query1 <- paste("local_res_atom_non_h_electron_sum ~ ", without_electron, sep = "")
electon_lm <- lm(joined_query1, numeric_data) #lm: fit linear models
electon_lm_summary <- summary(electon_lm)
Szacowana trafność regresji:
R^2
electon_lm_summary$r.squared
## [1] 0.0007736674
RMSE
sqrt(mean(electon_lm_summary$residuals^2))
## [1] 92.12367
Analogicznie postępujemy dla liczby elektronów:
without_atoms <- names(numeric_data)[-which(names(numeric_data) =='local_res_atom_non_h_count')]
joined_query1 <- paste("local_res_atom_non_h_count ~ ",without_atoms,sep = "")
atom_lm <- lm(joined_query1, numeric_data) #lm: fit linear models
atom_lm_summary <- summary(atom_lm)
Szacowana trafność regresji:
R^2
atom_lm_summary$r.squared
## [1] 0.0006942596
RMSE
sqrt(mean(atom_lm_summary$residuals^2))
## [1] 13.73096
Na początek oczyszczamy dane z kolumn, których nie powinniśmy brać do klasyfikacji:
columns_to_remove <- c('title','pdb_code','res_id',' chain_id','grid_space',' solvent_radius',' solvent_opening_radius',' part_step_FoFc_std_min',' part_step_FoFc_std_max',' part_step_FoFc_std_step')
cleaned_data <- data[ , -which(names(data) %in% columns_to_remove)]
cleaned_data <- cleaned_data %>% select(-starts_with("dict_"), -starts_with("local_"))
cleaned_data <- cleaned_data[, sapply(cleaned_data, nlevels) != 1]
Stratyfikujemy zbiór, trenujemy za pomocą algorytmu random forest:
Na koniec podsumowanie wyników:
# rfClasses <- predict(fit, newdata = testing)
# confusionMatrix(data = rfClasses, testing$Class)